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Abstract 

The Monte Carlo program CATCH (Capture And Transport of CHarged particles 
in a crystal) for the simulation of planar channelling in bent crystals is presented. The 
program tracks a charged particle through the distorted-crystal lattice with the use of 
continuous-potential approximation and the non-diffusion approach to the processes of 
scattering on electrons and nuclei. The output consists of the exit angular distributions, 
the energy loss spectra, and the spectra of any close-encounter process of interest. The 
curvature variability, face twist, and various surface imperfections of the real crystal 
can be taken into account. 
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1 Introduction 



Channelling of a beam of charged particles in a bent monocrystal, is going to become 
a working tool for the next generation of accelerators Therefore there is a need 
for a theoretical tool which describes a whole set of experimental data on channelling 
in the GeV range, and which also simulates the processes important for the future 
applications. 

The output should be the distribution of exiting primary and secondary particles, the 
energy loss in crystal, and any other interesting quantity related to channelling. Since 
these processes are sensitive to the orientation, the simulation should track every particle 
through the crystal lattice computing the probabihty of any process as a function of 
the co-ordinates. 

Historically, one of the independent discoveries of channelling, in the early sixties, 
was done in a Monte Carlo simulation of low-energy (<MeV) ions propagating in crys- 
tals 0. The very low ranges and the thin crystals used allowed the study of binary 
collisions of the incident ion with the atoms of the crystal. At GeV energy , crystals 
of a few centimetres in length are used, therefore tracking with binary collisions would 
take a considerable time. Instead, an approach with continuous potential introduced by 
J. Lindhard can be used. In this approach one considers collisions of the incoming 
particle with the atomic strings or planes instead of separate atoms, if the particle is 
sufficiently aligned with respect to the crystallographic axis or plane. The validity of 
doing this improves with the increase of the particle velocity 0. 

Besides the motion in the potential one must take into account the scattering. This 
makes it necessary to use either kinetic equations [Q, or computer simulation to 
transport particles through a crystal. The general feature of the methods described in 
Refs. 0, ^ is the use of the diffusion approach which omits the single scattering acts. 
However, in Monte Carlo methods it is easy to include the single collisions with nuclei 
and electrons. Besides their influence on channelling, these close-encounter processes are 
the source of secondary particles emitted from the crystal. Moreover, such collisions 
with electrons provide interesting peculiarities of the energy loss spectra in aligned 
crystals. Here we describe the Monte Carlo program CATCH for the simulation of 
planar channelling, which does not use the diffusion approximation. 

To simulate the nuclear collisions we use the LUND routines; the CATCH pro- 
gram serves as a frame to provide the orientational dependence of these processes. 
Another useful feature of the CATCH program is the various imperfections of the 
crystal surface, incorporated in the simulation because of their essential role in the 
crystal-assisted beam extraction. Version 1.4 of CATCH allows variable curvature of 
the crystal, both longitudinal and transversal, in order to set, for instance, a twist of 
crystal faces. 

2 Continuous potential 

For the potential of the atomic plane we use the Moliere approximation; details can be 
found e.g. in the review by D. Gemmel The static-lattice potential is modified to 
take into account the thermal vibrations of the lattice atoms; this is done by integration 
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over the Gaussian distribution of the atom displacement. Bending of the crystal has no 
effect on this potential. However, it causes a centrifugal force in the non-inertial frame 
related to the atomic planes. To solve the equation of motion in the potential U{x) of 
the bent crystal, as a first approximation to the transport of a particle, 

(fx dU{x) pv 
^'"d^ " d^ ~ ' ^ ' 

[x being the transversal, z the longitudinal coordinate, pv the particle longitudinal 
momentum and velocity product, R{z) the local radius of curvature), we use the fast 
form of the Verlet algorithm 

Xi+i -Xi = {di + Q.hfi8z)5z , (2) 

e,+^ -6^ = 0.5(/,+i + f,)6z (3) 

with 6 for dx/dz, f for the 'force', and 6z for the step. It was chosen over the other 
second order algorithms for non-linear equations of motion, such as Euler- Cromer's 
and Beeman's, owing to the better conservation of the transverse energy shown in the 
potential motion. 

Figure shows the simulated phase trajectories of protons on the plane {x, 9) in 
Si (111). Scattering was also included; as a result, there are no channelled protons 
near the atomic planes. With version 1.4 of CATCH one can specify in the input cards 
the planar geometry equivalent either to the Si(llO) geometry (equidistant planes), or 
to the Si(lll) one. It is easy to 'build' any other geometry by editing the program. 
The bending curvature and the crystal plane orientation can be arbitrary functions of 
spacial coordinates. Any data measured (or assumed) for the real crystal shape can be 
implemented in the simulation. 



3 Scattering 

Beam bending by a crystal is due to the trapping of some particles in the potential well 
[/(x), where they then follow the direction of the atomic planes. This simple picture is 
disturbed by scattering processes which could cause (as result of one or many acts) the 
trapped particle to come to a free state (feed out, or de channelling process), and an 
initially free particle to be trapped in the channelled state (feed in, or volume capture). 



3.1 Scattering on electrons 

Feed out is mostly due to scattering on electrons f^, because the channelled particles 
keep far from the nuclei. The mean energy loss in this scattering can be written as 
follows 0: 

dE D 2m,c^/?^^ 2 5 T^, 
~~dl^ ~I ^ ~ 2 ^-^P^^^l ' 

with D = An NAr'^rrieC^ z'^ ^p, z for the charge of the incident particle (in units of e), p 
for the crystal density, Z and A for atomic number and weight, and the other notation 
being standard 0. 
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The second part in the brackets is due to single colhsions and depends on the local 
density Pe{x) (normalized on the amorphous one) of electrons. The angle of scattering 
in soft acts can be computed as a random Gaussian with r.m.s. value 

^rms = -^{^E)soft, (5) 

where {5E)soft is the soft acts, contribution in Eq. (4). The probability of the hard 
collision contributing to the second part of Eq. (4) is computed at every step. The 
energy transfer T in such an act is generated according to the distribution function 
P{T): 

DnJx) 1 , , 

^(^) = 4^^- (6) 

The transverse momentum transfer q is equal to 



q=pm,T+iT/cf. (7) 
Its projections are used to modify the angles 9^ and 9y of the particle. 

3.2 Scattering on nuclei 

The multiple Coulomb scattering on nuclei is computed by the approximation Kitagawa- 
Ohtsuki [0: 

{^nud,scl = {9sc)amorph ' Pn{x) (8) 

i.e. the mean angle of scattering squared is proportional to the local density of nuclei 
p„(x). The density function is Gaussian with r.m.s. value u being the amplitude thermal 
vibration of the atom. 

The probability of nuclear collision, proportional to Pn{x), is checked at every step. 
If such a collision succeeds (then the flag ISTATU is set to -1), the LUND routine 
responsible for the event generation may be called. This scheme may be used for the 
description of any close-encounter process of interest; then one should define the collision 
length (ADLI in the input cards) properly. In version 1.4 of CATCH the secondary 
particles are normally not tracked through the crystal lattice, but their exit angles can 
be modified due to multiple scattering in the rest part of the crystal. In principle, a 
full-scale tracking of secondary particles (and even their products) is possible via the 
construction of the appropriate main program (description below) with loops and trees. 



4 Crystal imperfection 

The complete description of the program usage follows below. Here we show how the 
crystal imperfection is implemented. 

4.1 Variable curvature and twist 

If KRADV = is set in the input cards, the transverse curvature is set to 0. The 
longitudinal curvature is for the longitudinal coordinate z between and STRAIT 
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(length of the forthcoming straight part), and is constant equal to 1/RADlUS for z 
between STRAIT and STRAIT + DLINA (hence DLINA means the length of the bent 
part ). 

If KRADV = 1 is set in the input cards, the curvature at any point must be defined 
by the user in function CURVAT. The total length of the crystal is STRAIT + DLINA 
in any case. 

If KRADVY = is set in the input cards, the orientation of the crystal planes at 
both entry and exit faces is constant. If KRADVY = 1 is set in the input cards, these 
orientations are computed according to the functions TCFROY and TCENDY, for the 
front face and for the back face respectively. These functions must be supplied by the 
user. 

4.2 Surface effects 

In the case of beam extraction from the accelerator, extremely small impact parameters 
are possible, making the surface effects essential. The particle entering the crystal very 
close to its edge, can suffer from various additional factors: 

• miscut angle (between the atomic planes and the surface), 

• roughness (i.e. non flatness) of the surface, 

• possible amorphous layer, 

• bent surface. 

Therefore one must pay particular attention to the near-surface tracking, where the 
particule is entering and leaving the cristal materi al (due to the roughness, holes and 
bend), both coherent and non-coherent scattering in this peculiar region, bending in 
short channels, and so on. The surface effects mentioned above are simulated in the cur- 
rent version of CATCH. The roughness is expressed by a periodical function, asin(2;/A), 
where a is the amplitude of 'bumps' and A is their periodicity. The 'rough' crystalline 
material can be superimposed by a uniform amorphous layer. The position of surfaces 
is computed at every step in accordance with bending with the variable curvature. To 
set up these imperfections, simply give non-zero values to the corresponding parameters 
in the input cards. 

5 Usage of routines 

Figure 2 shows an example of usage of the CATCH routines. The subroutine XPREP 
asks for a filename with a problem description, and performs some preliminary pro- 
cedures. The subroutine XCATCH takes the entry values of ,x, Q^^ y, Qy and E and 
returns the exit values. The entry/exit values of X, Y, XP, YP and ENERGY are dou- 
ble precision . The units in CALL. ..(..) are radians and metres; the energy is returned 
modified according to the E loss in crystal; hitting/missing the crystal is checked, the 
flag for that is ISTATU (returned) : 
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• ISTATU = 0, crystal missed, 

• ISTATU = 1, or -1, crystal hit, 

• ISTATU < 0, there was nuclear interaction. 

In the case of interaction one may call the necessary routines to produce the secondary 
particles and to get their exit parameters modified for multiple scattering. The XPOST 
can be called to save the histograms, to fit the data, etc. Figure 3 shows an example of 
the input file. The input data is listed below with comments when necessary. 

• title 

• KGEOM, geometry type: either 110 (equidistant planes) or 111 (planes like those 
in the Diamond(lll) lattice) 

• KTRACE, trace level: 

— if < 0, no output during tracking 

— if = 1, some output at every step of tracking 

— if > 2, all current parameters are printed at every step 

• KRADV, switch on/off R = R{z) option; with KRADV = 1 one must sup- 
ply FUNCTION CURVAT, evaluating the curvature 1/R for every point; with 
KRADV = the curvature is 1/RADIUS. 

• KRADVY, switch on/off face twist; with KRADVY=1 one must supply functions 
TCFROY and TCENDY, evaluating the orientation of the planes for any point 
at the crystal faces; with KRADVY = that orientation is constant. 

• II, 12, initial numbers for pseudo-random generator RAN(I1,I2) 

• DZ, step size, microns. Try different DZ to be sure the result is independent of 
the step size. A step of the order of 1 fim is suggested for 450 GeV/c protons in 
Si(llO). The reasonable value of DZ should be scaled like VE (like the oscillation 
period in channelling) 

• ENERGY, GeV 

• PMASS, mass of the beam particle, (GeV) 

• RADIUS, in cm, for the bent part of the crystal. When KRADV = 1 is set, the 
RADIUS is used to define the histogram windows only. 

• DLINA, length of the bent part of the crystal (cm). The whole length of the 
crystal is STRAIT + DLINA. When KRADV = is set, the forthcoming part 
of the length STRAIT is assumed to be straight. The following part of the th 
DLINA is assumed to be bent with R = RADIUS. 

• STRAIT, length of the forthcoming straight part, (cm) 
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• HDL, interplanar half-spacing, for (110) or (lll)L, in A 

• HDS, interplanar half-spacing for (lll)S, in A. Must be for KGEOM=110 

• U, thermal vibration amplitude of the lattice atoms, (A) 

• DENS, crystal density, (gram/cm^) 

• RDLI, radia ion length, (cm) 

• ADLI, absorption length, (cm) 

• ZN, charge of the nucleus 

• AN, atomic weight 

• XOXTL, the crystal entry-face left-edge x-coordinate, (cm) 

• XIXTL, the crystal entry-face right-edge x-coordinate, (cm) 

• YOXTL, the crystal entry-face left-edge y-coordinate, (cm) 

• YIXTL, the crystal entry-face right-edge y-coordinate, (cm) 

• THXTAL, crystal angle 6*^, (//rad) 

• THYTAL, crystal angle Oy, (//rad) 

• SKINTH, amorphous skin thickness, (cm) 

Some amorphous skin with uniform thickness SKINTH can be superimposed over 
the surface. 

• ROUGH, (cm) the roughness (= non-fiatness) amplitude of every surface. It is 
described by some periodical function, therefore you should supply in the input 
both its amphtude ROUGH and period ROUGPD. It is assumed that material 
under the rough surface is crystalline. 

• ROUGPD, surface roughness period, (cm) 

• SKIMOZ, (^rad) There is the same skin, SKINTH + ROUGH, at the entry face. 
If you wish, you can also set some angular spread SKIMOZ of the crystal lattice 
in the entry- face skin (in this skin only!). The SKIMOZ means the sigma of 
Gaussian distribution. 

• DISCUT, miscut angle, (/xrad) Note that its sign is essential; when your impact 
is near the left edge of the crystal, the negative sign (DISCUT<0) is preferable. 

The position of all surfaces in space, as a function of 9^ and 9y misalignment of the 
crystal (THXTAL and THYTAL in the input cards), as a function of bending with 
variable radius, as well as a function of the surface 'rough structure', is computed at 
every step. 
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All variables read from the input cards are placed into COMMON blocks, and hence 
you can vary them from the main program during the execution (if you wish to scan 
or to optimize). Every particle (both channelled and non-channelled) is tracked in the 
crystal (X, Y are changed) and can leak out through any surface. Near the rough surface 
it is even possible to leak out and be caught again many times, i.e. the particle traverses 
sequentially the crystal bumps and vacuum (or amorphous skin) between. Hence, for 
particles touching any surface one has both Coulomb multiple and coherent collisions. 
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PROGRAM XMAIN 
C 

C 

c 

An EXAMPLE of the 

C Monte Carlo simulation 

C of the planar channelling in bent crystals 

C 

DOUBLE PRECISION X,Y,XP,YP, ENERGY 
C 

C . . . 

CALL XPREP ! preliminary procedures 

C . . . 

NSTAT = 100 

C statistics cycle up to NSTAT events 

C . . . 

DO 40 IS = 1, NSTAT 
C . . . 

XP = GAUS (0. ,20.E-6) 
X = 2.E-3+l.E-6*RAN(Il,I2) 
YP = 0.0 
Y = 0.0 
ENERGY =120. 
C . . . 

CALL XCATCH (X,XP,Y,YP,ISTATU, ENERGY) ! crystal simulation 
C . . . 

if (ISTATU.LT.O) TYPE interacted' 
C . . . 

TYPE *,IS, XP.YP, ISTATU 
C . . . 

40 CONTINUE ! end of cycle 

C . . . 

CALL XPOST ! some output 

C . . . 

C ========================================================= 

STOP 
END 



Figure 2: An example of usage of the CATCH routines 



9 



! pseudo-random start 
! -//- 
step size, microns 



7700., ! GeV 

.938, ! mass of particle. 



GeV 



LHC-EXTRACTION simulation by CATCH $KEYS 
KGEOM =110, ! geometry type 
KTRACE =0, ! trace level 

KRADV = 0, ! switch on/off R = R(x,y,z) option 
KRADVY =0, ! switch on/off face twist 

$END 

$OPTN 

11 = 84837, 

12 = 53463, 
DZ = 5 . , ! 

$END 
$BEAM 
ENERGY 
PMASS 
$END 
$CRYS 
RADIUS 
DLINA 
STRAIT 
HDL = 
HDS = 
U 

DENS = 
RDLI = 
ADLI = 
ZN = 
AN = 
$END 
$GEOM 
XOXTL =0.2, 
XIXTL =0.5, 
YOXTL = -.15, 
YIXTL = .15, 
THXTAL = .0, ! 
THYTAL = .0, ! 
SKINTH = O.E-5, 
ROUGH = O.E-4, 
ROUGPD = .01, 
SKIMOZ = 10. , 
DISCUT = 0. , 
$END 



14., ! 
28.09, 



l.e4, ! cm 

7 . , ! length of bent part , cm 
= 0., ! length of forthcoming straight part, cm 
.96, ! interplanar half-spacing, (110) or (lll)L 
.0, ! interplanar half -spacing, (lll)S, Angst. 
.075, ! thermal vibration amplitude, Angstroems 
2.33, ! density, gramm/cm3 
9.38, ! radiation length, cm 
45., ! absorbtion length, cm 
charge of nucleus 
! atomic weight 



min X of the Xtal-face (cm) 
max X of the Xtal-face (cm) 
min Y of the Xtal-face (cm) 
max Y of the Xtal-face (cm) 
xtal angle wrt X (microrad) 
xtal angle wrt Y (microrad) 
! amorphous skin thickness (cm) 
! surface roughness (cm) 
surface roughness period (cm) 
skin mozaicity (microrad) 
miscut angle (microrad) 



Figure 3: An example of the input file 



10 



